library(foreign)
library(MASS)
library(pscl)
library(DataCombine)
library(stargazer)
library(plm)
library(lfe)
library(survival)
library(doBy)
library(plyr)

setwd("~/My Drive/Atrocity development/Conditional acceptance_10_19_21/Replication/")
main.grid <- read.dta("grid_atrocity_10_7_21.dta")
summary(main.grid)

#Number of countries
length(unique(na.omit(main.grid$gwno)))
#Number of cells per year
main.grid$one <- 1
cells <- summaryBy(one ~ gid, FUN=c("sum"), data=main.grid, keep.names = T)
#Divide by country
length(cells$one)/length(unique(na.omit(main.grid$gwno)))
###################
###################
###Main Analyses###
###################
###################

#########
###OLS###
#########
###Baseline
fel1 <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
               laglogpop|year+gid|0|gid, data=main.grid)

summary(fel1)
###Medium
fel2 <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
               laglogpop+excluded+p_polity2|year+gid|0|gid, data=main.grid)

summary(fel2)
###Full
fel3 <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
               laglogpop+lagnlights_mean+
               excluded+p_polity2 +
               logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3)
###Export to LaTex
stargazer(fel1,fel2,fel3)


#########
###GMM###
#########

####Baseline
##Keep only variables of interest
omit.base <- na.omit(main.grid[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid")])
#N units
N1 <- summaryBy(nlights_mean~gid, FUN=c("sum"), data=omit.base, keep.names = T)
length(N1$nlights_mean)
#Convert data to pgmm formal
pgmm.base <- pdata.frame(omit.base, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.1 <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                laglogpop| 
                lag(lagall_perp_ml, 2:7),
              data = pgmm.base, 
              effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.1)

####Medium
##Keep only variables of interest
omit.med <- na.omit(main.grid[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","excluded","p_polity2")])
#N units
N2 <- summaryBy(nlights_mean~gid, FUN=c("sum"), data=omit.med, keep.names = T)
length(N2$nlights_mean)
#Convert data to pgmm formal
pgmm.med <- pdata.frame(omit.med, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.2 <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                laglogpop+excluded+p_polity2| 
                lag(lagall_perp_ml, 2:7),
              data = pgmm.med, 
              effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.2)


####Full
##Keep only variables of interest
omit.full <- na.omit(main.grid[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#N units
N3 <- summaryBy(nlights_mean~gid, FUN=c("sum"), data=omit.full, keep.names = T)
length(N3$nlights_mean)
#Convert data to pgmm formal
pgmm.full <- pdata.frame(omit.full, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3 <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                laglogpop+lagnlights_mean+
                excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                lag(lagall_perp_ml, 2:7),
              data = pgmm.full, 
              effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3)


###Export to Latex
stargazer(gmm.1, gmm.2, gmm.3)




###Export all to Latex
stargazer(fel1, fel2, fel3, gmm.1, gmm.2, gmm.3)


###############################
###############################
###Size of Chnage Thresholds###
###############################
###############################

###Create threshold indicators
#One year lag
main.grid$lagall_perp_ml_one <- ifelse(main.grid$lagall_perp_ml>0,1,0)
table(main.grid$lagall_perp_ml_one)
main.grid$lagall_perp_ml_ten <- ifelse(main.grid$lagall_perp_ml>9,1,0)
table(main.grid$lagall_perp_ml_ten)
main.grid$lagall_perp_ml_twenty <- ifelse(main.grid$lagall_perp_ml>19,1,0)
table(main.grid$lagall_perp_ml_twenty)
#Two year lag
main.grid$lag2all_perp_ml_one <- ifelse(main.grid$lag2all_perp_ml>0,1,0)
table(main.grid$lag2all_perp_ml_one)
main.grid$lag2all_perp_ml_ten <- ifelse(main.grid$lag2all_perp_ml>9,1,0)
table(main.grid$lag2all_perp_ml_ten)
main.grid$lag2all_perp_ml_twenty <- ifelse(main.grid$lag2all_perp_ml>19,1,0)
table(main.grid$lag2all_perp_ml_twenty)


################
###Run Models###
################
###One
fel3.one <- felm(nlights_mean~lagall_perp_ml_one+lag2all_perp_ml_one+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +
                   logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.one)
##Extract relevant quantities
#One year lag
mean.1 <- (summary(fel3.one)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.1  <- (summary(fel3.one)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.one)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.1 <- (summary(fel3.one)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) -1.97*((summary(fel3.one)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#Two year lag
mean.1b <- (summary(fel3.one)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.1b  <- (summary(fel3.one)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.one)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.1b <- (summary(fel3.one)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) -1.97*((summary(fel3.one)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))


###Ten
fel3.ten <- felm(nlights_mean~lagall_perp_ml_ten+lag2all_perp_ml_ten+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +
                   logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.ten)
##Extract relevant quantities
#One year lag
mean.10 <- (summary(fel3.ten)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.10  <- (summary(fel3.ten)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.ten)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.10 <- (summary(fel3.ten)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.ten)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#Two year lag
mean.10b <- (summary(fel3.ten)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.10b  <- (summary(fel3.ten)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.ten)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.10b <- (summary(fel3.ten)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.ten)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))


###Twenty
fel3.twenty <- felm(nlights_mean~lagall_perp_ml_twenty+lag2all_perp_ml_twenty+
                      laglogpop+lagnlights_mean+
                      excluded+p_polity2 +
                      logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.twenty)
##Extract relevant quantities
#One year lag
mean.20 <- (summary(fel3.twenty)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.20  <- (summary(fel3.twenty)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.twenty)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.20 <- (summary(fel3.twenty)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.twenty)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#TWo year lag
mean.20b <- (summary(fel3.twenty)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.20b  <- (summary(fel3.twenty)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.twenty)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.20b <- (summary(fel3.twenty)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.twenty)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))


#
stargazer(fel3.one,fel3.ten,fel3.twenty)

#Ntl average (without zeros)
main.grid$nlights_mean_na <-  ifelse(main.grid$nlights_mean>0, (main.grid$nlights_mean), NA)
summary(main.grid$nlights_mean_na)
#Median
median(main.grid$nlights_mean_na, na.rm=T)
#Mean
mean(main.grid$nlights_mean_na, na.rm=T)


################
###Run Models###
################
###One
fel3.one <- felm(nlights_mean~lagall_perp_ml_one+lag2all_perp_ml_one+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +
                   logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.one)
##Extract relevant quantities
#One year lag
mean.1 <- (summary(fel3.one)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.1  <- (summary(fel3.one)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.one)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.1 <- (summary(fel3.one)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) -1.97*((summary(fel3.one)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#Two year lag
mean.1b <- (summary(fel3.one)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.1b  <- (summary(fel3.one)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.one)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.1b <- (summary(fel3.one)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) -1.97*((summary(fel3.one)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))


###Ten
fel3.ten <- felm(nlights_mean~lagall_perp_ml_ten+lag2all_perp_ml_ten+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +
                   logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.ten)
##Extract relevant quantities
#One year lag
mean.10 <- (summary(fel3.ten)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.10  <- (summary(fel3.ten)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.ten)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.10 <- (summary(fel3.ten)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.ten)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#Two year lag
mean.10b <- (summary(fel3.ten)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.10b  <- (summary(fel3.ten)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.ten)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.10b <- (summary(fel3.ten)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.ten)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))


###Twenty
fel3.twenty <- felm(nlights_mean~lagall_perp_ml_twenty+lag2all_perp_ml_twenty+
                      laglogpop+lagnlights_mean+
                      excluded+p_polity2 +
                      logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.twenty)
##Extract relevant quantities
#One year lag
mean.20 <- (summary(fel3.twenty)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.20  <- (summary(fel3.twenty)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.twenty)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.20 <- (summary(fel3.twenty)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.twenty)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#TWo year lag
mean.20b <- (summary(fel3.twenty)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.20b  <- (summary(fel3.twenty)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.twenty)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.20b <- (summary(fel3.twenty)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.twenty)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))


#
stargazer(fel3.one,fel3.ten,fel3.twenty)

#Ntl average (without zeros)
main.grid$nlights_mean_na <-  ifelse(main.grid$nlights_mean>0, (main.grid$nlights_mean), NA)
summary(main.grid$nlights_mean_na)
#Median
median(main.grid$nlights_mean_na, na.rm=T)
#Mean
mean(main.grid$nlights_mean_na, na.rm=T)

#########################
###Plot Change in Odds###
#########################

##################
###ONE YEAR LAG###

###Set up plotting space
par(mfrow=c(1,1))
par(mar=c(5,4,2,2) + 0.1)
par(cex.lab=1)
par(cex.axis=1)
par(cex.main=1)

#Set up axes
ylim<-c(-200,30)
xlim<-c(-.1,21)
x<-c(1,10,20)
y<-c(mean.1*100,mean.10*100,mean.20*100)
yupper<-c(upper.1*100,upper.10*100,upper.20*100)
ylower<-c(lower.1*100,lower.10*100,lower.20*100)

#Create estimates for a trendline
fit <- lm(y~x)

###plot actual effects
setwd("~/My Drive/Atrocity development/")
#open file
jpeg("changethresh1.jpeg", width = 6, height = 6, units = 'in', res = 500)
#plot values
plot(x,y,ylim=ylim,xlim=xlim,col="Red",pch=16,lwd=2, 
     axes=FALSE, ann=FALSE,cex=1.5)

#add 95% CI's
segments(x,yupper,x,ylower,col="Red",lwd=2,cex=1.5)

#Add horizontal line to plot
abline(h=0,lty=3)

#Add a trendline
lines(x, fitted(fit), col="blue")

#Label X-axis values
axis(1, at=c(1,10,20), lab=c("1 atrocity","10 atrocities","20 atrocities"), tick = TRUE)

#Add a title
title(ylab="% Decrease in socioeconomic activity (from sample mean)")

#Adjust Y-Axis Labels
axis(2, at=c((-200), (-100), 0), lab=c("-200%","-100%","0%"))

#Add box around plot
box()
#close file
dev.off()

##################
###TWO YEAR LAG###

###Set up plotting space
par(mfrow=c(1,1))
par(mar=c(5,4,2,2) + 0.1)
par(cex.lab=1)
par(cex.axis=1)
par(cex.main=1)

#Set up axes
ylim<-c(-270,5)
xlim<-c(-.1,21)
x<-c(1,10,20)
y<-c(mean.1b*100,mean.10b*100,mean.20b*100)
yupper<-c(upper.1b*100,upper.10b*100,upper.20b*100)
ylower<-c(lower.1b*100,lower.10b*100,lower.20b*100)

#Create estimates for a trendline
fit <- lm(y~x)

###plot actual effects
setwd("~/My Drive/Atrocity development/")
#open file
jpeg("changethresh2.jpeg", width = 6, height = 6, units = 'in', res = 500)
#plot values
plot(x,y,ylim=ylim,xlim=xlim,col="Red",pch=16,lwd=2, 
     axes=FALSE, ann=FALSE,cex=1.5)

#add 95% CI's
segments(x,yupper,x,ylower,col="Red",lwd=2,cex=1.5)

#Add horizontal line to plot
abline(h=0,lty=3)

#Add a trendline
lines(x, fitted(fit), col="blue")

#Label X-axis values
axis(1, at=c(1,10,20), lab=c("1 atrocity","10 atrocities","20 atrocities"), tick = TRUE)

#Add a title
title(ylab="% Decrease in socioeconomic activity (from sample mean)")

#Adjust Y-Axis Labels
axis(2, at=c((-300), (-200), (-100), 0), lab=c("-300%","-200%","-100%","0%"))

#Add box around plot
box()
#close file
dev.off()


##############################
##############################
###Non Lags and Deeper Lags###
##############################
##############################

################
###OLS Models###
################

###Lag 0-1
fel3.lag0_1 <- felm(nlights_mean~all_perp_ml+lagall_perp_ml+
                      laglogpop+lagnlights_mean+
                      excluded+p_polity2 +
                      logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.lag0_1)

###Lag 0-2
fel3.lag0_2 <- felm(nlights_mean~all_perp_ml+lagall_perp_ml+lag2all_perp_ml+
                      laglogpop+lagnlights_mean+
                      excluded+p_polity2 +
                      logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.lag0_2)


###3rd lag
fel3.lag3 <- felm(nlights_mean~lag3all_perp_ml+
                    laglogpop+lagnlights_mean+
                    excluded+p_polity2 +
                    logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.lag3)

##Fourth lag
fel3.lag4 <- felm(nlights_mean~lag(all_perp_ml,4)+
                    laglogpop+lagnlights_mean+
                    excluded+p_polity2 +
                    logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.lag4)

###Lag 1-3
fel3.lag1_3 <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+lag3all_perp_ml+
                      laglogpop+lagnlights_mean+
                      excluded+p_polity2 +
                      logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.lag1_3)
###Lag 1-4
fel3.lag1_4 <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+lag3all_perp_ml+lag(all_perp_ml,4)+
                      laglogpop+lagnlights_mean+
                      excluded+p_polity2 +
                      logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.lag1_4)

###Lag 0-3
fel3.lag0_3 <- felm(nlights_mean~all_perp_ml+lagall_perp_ml+lag2all_perp_ml+lag3all_perp_ml+
                      laglogpop+lagnlights_mean+
                      excluded+p_polity2 +
                      logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.lag0_3)


##OLS
#stargazer(fel3.lag0_1,fel3.lag0_2,fel3.lag3,fel3.lag4,fel3.lag1_3,fel3.lag1_4,fel3.lag0_3)


################
###OLS Models###
################


###Lag 0-1
##Keep only variables of interest
omit.full.lag0_1 <- na.omit(main.grid[,c("nlights_mean","all_perp_ml","lagall_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#Convert data to pgmm formal
pgmm.full.lag0_1 <- pdata.frame(omit.full.lag0_1, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.lag0_1 <- pgmm(nlights_mean~all_perp_ml+lagall_perp_ml+
                       laglogpop+lagnlights_mean+
                       excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                       lag(all_perp_ml, 2:7),
                     data = pgmm.full.lag0_1, 
                     effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.lag0_1)


###Lag 0-2
##Keep only variables of interest
omit.full.lag0_2 <- na.omit(main.grid[,c("nlights_mean","all_perp_ml","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#Convert data to pgmm formal
pgmm.full.lag0_2 <- pdata.frame(omit.full.lag0_2, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.lag0_2 <- pgmm(nlights_mean~all_perp_ml+lagall_perp_ml+lag2all_perp_ml+
                       laglogpop+lagnlights_mean+
                       excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                       lag(lagall_perp_ml, 2:7),
                     data = pgmm.full.lag0_2, 
                     effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.lag0_2)

###3rd lag
##Keep only variables of interest
omit.full.lag3 <- na.omit(main.grid[,c("nlights_mean","lag3all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#Convert data to pgmm formal
pgmm.full.lag3 <- pdata.frame(omit.full.lag3, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.lag3 <- pgmm(nlights_mean~lag3all_perp_ml+
                     laglogpop+lagnlights_mean+
                     excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                     lag(lag3all_perp_ml, 2:7),
                   data = pgmm.full.lag3, 
                   effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.lag3)

##Fourth lag
###Run model
gmm.3.lag4 <- pgmm(nlights_mean~lag(lag3all_perp_ml,1)+
                     laglogpop+lagnlights_mean+
                     excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                     lag(lag3all_perp_ml, 3:8),
                   data = pgmm.full.lag3, 
                   effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.lag4)

###Lag 1-3
##Keep only variables of interest
omit.full.lag1_3 <- na.omit(main.grid[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","lag3all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#Convert data to pgmm formal
pgmm.full.lag1_3 <- pdata.frame(omit.full.lag1_3, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.lag1_3 <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+lag3all_perp_ml+
                       laglogpop+lagnlights_mean+
                       excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                       lag(lag3all_perp_ml, 2:7),
                     data = pgmm.full.lag1_3, 
                     effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.lag1_3)


###Lag 1-4
###Run model
gmm.3.lag1_4 <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+lag3all_perp_ml+lag(lag3all_perp_ml,1)+
                       laglogpop+lagnlights_mean+
                       excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                       lag(lag3all_perp_ml, 3:8),
                     data = pgmm.full.lag1_3, 
                     effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.lag1_4)


###Lag 0-3
##Keep only variables of interest
omit.full.lag0_3 <- na.omit(main.grid[,c("nlights_mean","all_perp_ml","lagall_perp_ml","lag2all_perp_ml","lag3all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#Convert data to pgmm formal
pgmm.full.lag0_3 <- pdata.frame(omit.full.lag0_3, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.lag0_3 <- pgmm(nlights_mean~all_perp_ml+lagall_perp_ml+lag2all_perp_ml+lag3all_perp_ml+
                       laglogpop+lagnlights_mean+
                       excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                       lag(lag3all_perp_ml, 2:7),
                     data = pgmm.full.lag0_3, 
                     effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.lag0_3)


###Export both OLS and GMMs to LaTex
stargazer(fel3.lag0_1,fel3.lag0_2,fel3.lag3,fel3.lag4,fel3.lag1_3,fel3.lag1_4,fel3.lag0_3,gmm.3.lag0_1,gmm.3.lag0_2,gmm.3.lag3,gmm.3.lag4,gmm.3.lag1_3,gmm.3.lag1_4,gmm.3.lag0_3)


###########################################################
###########################################################
###Plot Mean Coefficients To Determine Over Time Effects###
###########################################################
###########################################################


######################
###OLS Coefficients###
######################

##Average the coefficients
at0 <- ((0.005)+(0.009)+0.012)/3
at1 <- ((-0.040)+(-0.032)+(-0.034)+(-0.036)+(-0.036))/5
at2 <- ((-0.038)+(-0.045)+(-0.046)+(-0.046))/4
at3 <- ((0.009)+(0.020)+(0.019)+(0.019))/4
at4 <- ((-0.005)+(0.012))/2

#Create data frame 
at.comb <- data.frame(rbind(at0,at1,at2,at3,at4))
at.comb$lag <- c("0","1","2","3","4")
names(at.comb) <- c("Mean Coeff.", "Lag")
#Normalize by percent effect
at.comb$coef_norm <- (at.comb$`Mean Coeff.`/(1.192))*100

###Set up plotting space
par(mfrow=c(1,1))
par(mar=c(5,4,2,2) + 0.1)
par(cex.lab=1)
par(cex.axis=1)
par(cex.main=1)

#Set up axes
ylim<-c(-5,2.5)
xlim<-c(-.1,4.1)
x<-c(0,1,2,3,4)
y<-at.comb$coef_norm

#Create estimates for a trendline
fit <- lm(y~x)

###plot actual effects
setwd("~/My Drive/Atrocity development/")
#open file
jpeg("olseff.jpeg", width = 6, height = 6, units = 'in', res = 500)
#plot values
plot(x,y,ylim=ylim,xlim=xlim,col="Red",pch=16,lwd=2, 
     axes=FALSE, ann=FALSE,cex=1.5)

#add 95% CI's
#segments(x,yupper,x,ylower,col="Red",lwd=2,cex=1.5)

#Add horizontal line to plot
abline(h=0,lty=3)

#Add a trendline
lines(x, fitted(fit), col="blue")

#Label X-axis values
axis(1, at=c(0,1,2,3,4), lab=c("No lag", "1 year","2 years","3 years", "4 years"), tick = TRUE)

#Add a title
title(ylab="Average Model Effects (%)")

#Adjust Y-Axis Labels
axis(2, at=c((-5), (-2.5), 0, 2.5), lab=c("-5%","-2.5%","0%", "2.5%"))

#Add box around plot
box()
#close file
dev.off()


######################
###GMM Coefficients###
######################

##Average the coefficients
at0 <- ((-0.058)+(-0.028)+0.039)/3
at1 <- ((-0.057)+(-0.112)+(-0.073)+(-0.165)+(-0.041))/5
at2 <- ((-0.109)+(-0.159)+(-0.290)+(-0.148))/4
at3 <- ((-0.275)+(-0.260)+(-0.552)+(-0.258))/4
at4 <- ((-0.185)+(0.44))/2

#Create data frame 
at.comb <- data.frame(rbind(at0,at1,at2,at3,at4))
at.comb$lag <- c("0","1","2","3","4")
names(at.comb) <- c("Mean Coeff.", "Lag")
#Normalize by percent effect
at.comb$coef_norm <- (at.comb$`Mean Coeff.`/(1.192))*100

###Set up plotting space
par(mfrow=c(1,1))
par(mar=c(5,4,2,2) + 0.1)
par(cex.lab=1)
par(cex.axis=1)
par(cex.main=1)

#Set up axes
ylim<-c(-40,20)
xlim<-c(-.1,4.1)
x<-c(0,1,2,3,4)
y<-at.comb$coef_norm

#Create estimates for a trendline
fit <- lm(y~x)

###plot actual effects
setwd("~/My Drive/Atrocity development/")
#open file
jpeg("gmmeff.jpeg", width = 6, height = 6, units = 'in', res = 500)
#plot values
plot(x,y,ylim=ylim,xlim=xlim,col="Red",pch=16,lwd=2, 
     axes=FALSE, ann=FALSE,cex=1.5)

#add 95% CI's
#segments(x,yupper,x,ylower,col="Red",lwd=2,cex=1.5)

#Add horizontal line to plot
abline(h=0,lty=3)

#Add a trendline
lines(x, fitted(fit), col="blue")

#Label X-axis values
axis(1, at=c(0,1,2,3,4), lab=c("No lag", "1 year","2 years","3 years", "4 years"), tick = TRUE)

#Add a title
title(ylab="Average Model Effects (%)")

#Adjust Y-Axis Labels
axis(2, at=c((-40), (-20), 0, 20), lab=c("-40%","-20%","0%", "20%"))

#Add box around plot
box()
#close file
dev.off()


###########################
###########################
###High Importance Areas###
###########################
###########################

###########
###Urban###
###########
#Subset relevant data
main.grid.urb <- subset(main.grid, urban_gc>0)

###OLS
fel3.urb <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +
                   logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid.urb)
summary(fel3.urb)
#One year lag
mean.urb <- (summary(fel3.urb)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.urb  <- (summary(fel3.urb)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.urb)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.urb <- (summary(fel3.urb)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.urb)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#Two year lag
mean.urbb <- (summary(fel3.urb)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.urbb  <- (summary(fel3.urb)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.urb)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.urbb <- (summary(fel3.urb)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.urb)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))


###GMM
##Keep only variables of interest
omit.full.urb <- na.omit(main.grid.urb[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#omit.full$c_gid <- cluster(omit.full$gid)
#Convert data to pgmm formal
pgmm.full.urb <- pdata.frame(omit.full.urb, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.urb <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                    laglogpop+lagnlights_mean+
                    excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                    lag(lagall_perp_ml, 2:7),
                  data = pgmm.full.urb, 
                  effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.urb)



##############
###Big City###
##############
#Subset data
main.grid.city <- subset(main.grid, ttime_mean<=120)

###OLS
fel3.city <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                    laglogpop+lagnlights_mean+
                    excluded+p_polity2 +
                    logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid.city)

summary(fel3.city)
#One year lag
mean.city <- (summary(fel3.city)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.city  <- (summary(fel3.city)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.city)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.city <- (summary(fel3.city)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.city)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#Two year lag
mean.cityb <- (summary(fel3.city)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.cityb  <- (summary(fel3.city)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.city)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.cityb <- (summary(fel3.city)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.city)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))



###GMM
##Keep only variables of interest
omit.full.city <- na.omit(main.grid.city[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#omit.full$c_gid <- cluster(omit.full$gid)
#Convert data to pgmm formal
pgmm.full.city <- pdata.frame(omit.full.city, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.city <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                     laglogpop+lagnlights_mean+
                     excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                     lag(lagall_perp_ml, 2:7),
                   data = pgmm.full.city, 
                   effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.city)
#One year lag
mean.city <- (summary(fel3.city)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.city  <- (summary(fel3.city)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.city)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.city <- (summary(fel3.city)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.city)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#Two year lag
mean.cityb <- (summary(fel3.city)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.cityb  <- (summary(fel3.city)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.city)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.cityb <- (summary(fel3.city)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.city)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))



##################
###Capital City###
##################
#Subset relevant data
main.grid.cap <- subset(main.grid, capdist<=100)

###OLS
fel3.cap <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +
                   logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid.cap)

summary(fel3.cap)
#One year lag
mean.cap <- (summary(fel3.cap)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.cap  <- (summary(fel3.cap)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.cap)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.cap <- (summary(fel3.cap)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.cap)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#Two year lag
mean.capb <- (summary(fel3.cap)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.capb  <- (summary(fel3.cap)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3.cap)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.capb <- (summary(fel3.cap)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3.cap)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))


###GMM
##Keep only variables of interest
omit.full.cap <- na.omit(main.grid.cap[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#omit.full$c_gid <- cluster(omit.full$gid)
#Convert data to pgmm formal
pgmm.full.cap <- pdata.frame(omit.full.cap, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.cap <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                    laglogpop+lagnlights_mean+
                    excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                    lag(lagall_perp_ml, 2:7),
                  data = pgmm.full.cap, 
                  effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.cap)

###Export both OLS and GMMs to LaTex
stargazer(fel3.urb,gmm.3.urb,fel3.city,gmm.3.city,fel3.cap,gmm.3.cap)


###Create measures for the regular model
#One year lag
mean.reg <- (summary(fel3)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T)
upper.reg  <- (summary(fel3)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.reg <- (summary(fel3)$coefficients[1])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3)$coefficients[1,2])/mean(main.grid$nlights_mean, na.rm=T))
#Two year lag
mean.regb <- (summary(fel3)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T)
upper.regb  <- (summary(fel3)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) + 1.97*((summary(fel3)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))
lower.regb <- (summary(fel3)$coefficients[2])/mean(main.grid$nlights_mean, na.rm=T) - 1.97*((summary(fel3)$coefficients[2,2])/mean(main.grid$nlights_mean, na.rm=T))




#########################
###Plot Change in Odds###
#########################

##################
###ONE YEAR LAG###

###Set up plotting space
par(mfrow=c(1,1))
par(mar=c(5,4,2,2) + 0.1)
par(cex.lab=1)
par(cex.axis=1)
par(cex.main=1)

#Set up axes
ylim<-c(-7,1)
xlim<-c(-.1,31)
x<-c(1,10,20,30)
y<-c(mean.reg*100,mean.urb*100,mean.city*100,mean.cap*100)
yupper<-c(upper.reg*100,upper.urb*100,upper.city*100,upper.cap*100)
ylower<-c(lower.reg*100,lower.urb*100,lower.city*100,lower.cap*100)

#Create estimates for a trendline
fit <- lm(y~x)

###plot actual effects
setwd("~/My Drive/Atrocity development/")
#open file
jpeg("changethreshUrb.jpeg", width = 6, height = 6, units = 'in', res = 500)
#plot values
plot(x,y,ylim=ylim,xlim=xlim,col="Red",pch=16,lwd=2, 
     axes=FALSE, ann=FALSE,cex=1.5)

#add 95% CI's
segments(x,yupper,x,ylower,col="Red",lwd=2,cex=1.5)

#Add horizontal line to plot
abline(h=0,lty=3)

#Add a trendline
lines(x, fitted(fit), col="blue")

#Label X-axis values
axis(1, at=c(1,10,20,30), lab=c("All", "Urban","City","Capital"), tick = TRUE)

#Add a title
title(ylab="% Decrease in socioeconomic activity (from sample mean)")

#Adjust Y-Axis Labels
axis(2, at=c((-7), (-3.5), 0), lab=c("-7%","-3.5%","0%"))

#Add box around plot
box()
#close file
dev.off()



##################
###TWO YEAR LAG###

###Set up plotting space
par(mfrow=c(1,1))
par(mar=c(5,4,2,2) + 0.1)
par(cex.lab=1)
par(cex.axis=1)
par(cex.main=1)

#Set up axes
ylim<-c(-7,1)
xlim<-c(-.1,31)
x<-c(1,10,20,30)
y<-c(mean.regb*100,mean.urbb*100,mean.cityb*100,mean.capb*100)
yupper<-c(upper.regb*100,upper.urbb*100,upper.cityb*100,upper.capb*100)
ylower<-c(lower.regb*100,lower.urbb*100,lower.cityb*100,lower.capb*100)

#Create estimates for a trendline
fit <- lm(y~x)

###plot actual effects
setwd("~/My Drive/Atrocity development/")
#open file
jpeg("changethreshUrb2.jpeg", width = 6, height = 6, units = 'in', res = 500)
#plot values
plot(x,y,ylim=ylim,xlim=xlim,col="Red",pch=16,lwd=2, 
     axes=FALSE, ann=FALSE,cex=1.5)

#add 95% CI's
segments(x,yupper,x,ylower,col="Red",lwd=2,cex=1.5)

#Add horizontal line to plot
abline(h=0,lty=3)

#Add a trendline
lines(x, fitted(fit), col="blue")

#Label X-axis values
axis(1, at=c(1,10,20,30), lab=c("All", "Urban","City","Capital"), tick = TRUE)

#Add a title
title(ylab="% Decrease in socioeconomic activity (from sample mean)")

#Adjust Y-Axis Labels
axis(2, at=c((-7), (-3.5), 0), lab=c("-7%","-3.5%","0%"))

#Add box around plot
box()
#close file
dev.off()



#######################
#######################
###Robustness Models###
#######################
#######################

####################
###Added Controls###
####################
###OLS
fel3.cont <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                    laglogpop+lagnlights_mean+
                    excluded+p_polity2 +
                    logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod+
                    droughtyr_spi+logprec_gpcp+temp+ucdp_type3_bin+lagucdp_type3_bin|year+gid|0|gid, data=main.grid)

summary(fel3.cont)
###GMM
#Keep only variables of interest
omit.full.2 <- na.omit(main.grid[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2","droughtyr_spi","logprec_gpcp","temp","ucdp_type3_bin","lagucdp_type3_bin")])
#Convert data to pgmm formal
pgmm.full.2 <- pdata.frame(omit.full.2, index=c("gid", "year"))

#set seed
set.seed(50)
#Run model
gmm.3.cont <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                     laglogpop+lagnlights_mean+
                     excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod+
                     droughtyr_spi+logprec_gpcp+temp+ucdp_type3_bin+lagucdp_type3_bin| 
                     lag(lagall_perp_ml, 2:7),
                   data = pgmm.full.2, 
                   effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.cont)


#######################
###Extended Controls###
#######################
#Create a log version of arm imports
main.grid$logwdi_armimp <- log(main.grid$wdi_armimp+1)
###OLS
fel3.cont.2 <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                      laglogpop+lagnlights_mean+
                      excluded+p_polity2 +
                      logwdi_expmil + logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod+
                      droughtyr_spi+logprec_gpcp+temp+ucdp_type3_bin+lagucdp_type3_bin+
                      logcapdist+logwdi_armimp+wdi_lifexp+ciri_empinx_new|year+gid|0|gid, data=main.grid)

summary(fel3.cont.2)

###GMM
#Keep only variables of interest
omit.full.3 <- na.omit(main.grid[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2","droughtyr_spi","logprec_gpcp","temp","ucdp_type3_bin","lagucdp_type3_bin","logwdi_expmil","logwdi_armimp","wdi_lifexp","ciri_empinx_new")])
#Convert data to pgmm formal
pgmm.full.3 <- pdata.frame(omit.full.3, index=c("gid", "year"))

#set seed
set.seed(50)
#Run model
gmm.3.cont.2 <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                       laglogpop+lagnlights_mean+
                       excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod+
                       droughtyr_spi+logprec_gpcp+temp+ucdp_type3_bin+lagucdp_type3_bin+
                       logwdi_armimp+wdi_lifexp+ciri_empinx_new| 
                       lag(lagall_perp_ml, 2:7),
                     data = pgmm.full.3, 
                     effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.cont.2)

####################
###Only NTL areas###
####################
#Subset sample to only places with NTL
main.grid.ntl <- subset(main.grid, nlights_mean>0)

###OLS
fel3.ntl <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +
                   logwdi_expmil + logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid.ntl)

summary(fel3.ntl)


###GMM
#Keep only variables of interest
omit.full.ntl <- na.omit(main.grid.ntl[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#Convert data to pgmm formal
pgmm.full.ntl <- pdata.frame(omit.full.ntl, index=c("gid", "year"))

#set seed
set.seed(50)

#Run model
gmm.3.ntl <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                    laglogpop+lagnlights_mean+
                    excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                    lag(lagall_perp_ml, 2:7),
                  data = pgmm.full.ntl, 
                  effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.ntl)


##########################
###Only Civil War Areas###
##########################
#Subset sample to only places with NTL
main.grid.cw <- subset(main.grid, lagucdp_type3_bin>0)

###OLS
fel3.cw <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                  laglogpop+lagnlights_mean+
                  excluded+p_polity2 +
                  logwdi_expmil + logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid.cw)

summary(fel3.cw)

###GMM
#Keep only variables of interest
omit.full.cw <- na.omit(main.grid.cw[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#Convert data to pgmm formal
pgmm.full.cw <- pdata.frame(omit.full.cw, index=c("gid", "year"))

#set seed
set.seed(50)

#Run model
gmm.3.cw <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                   lag(lagall_perp_ml, 2:7),
                 data = pgmm.full.cw, 
                 effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.cw)

############################################
###Only Atrocities by Formal State Actors###
############################################
###OLS
fel3.s <- felm(nlights_mean~lags_perp_ml+lag2s_perp_ml+
                 laglogpop+lagnlights_mean+
                 excluded+p_polity2 +
                 logwdi_expmil + logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.s)

###GMM
#Keep only variables of interest
omit.full.s <- na.omit(main.grid[,c("nlights_mean","lags_perp_ml","lag2s_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#Convert data to pgmm formal
pgmm.full.s <- pdata.frame(omit.full.s, index=c("gid", "year"))

#set seed
set.seed(50)

#Run model
gmm.3.s <- pgmm(nlights_mean~lags_perp_ml+lag2s_perp_ml+
                  laglogpop+lagnlights_mean+
                  excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                  lag(lags_perp_ml, 2:7),
                data = pgmm.full.s, 
                effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.s)



#########################################
###Only Atrocities by Non-State Actors###
#########################################
###OLS
fel3.ns <- felm(nlights_mean~lagns_perp_ml+lag2ns_perp_ml+
                  laglogpop+lagnlights_mean+
                  excluded+p_polity2 +
                  logwdi_expmil + logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.ns)


###GMM
#Keep only variables of interest
omit.full.ns <- na.omit(main.grid[,c("nlights_mean","lagns_perp_ml","lag2ns_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#Convert data to pgmm formal
pgmm.full.ns <- pdata.frame(omit.full.ns, index=c("gid", "year"))

#set seed
set.seed(50)

#Run model
gmm.3.ns <- pgmm(nlights_mean~lagns_perp_ml+lag2ns_perp_ml+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                   lag(lagns_perp_ml, 2:7),
                 data = pgmm.full.ns, 
                 effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.ns)


###Export to LaTex
##OLS
stargazer(fel3.cont,fel3.cont.2,fel3.ntl,fel3.cw,fel3.orig,fel3.s,fel3.ns)
##GMM
stargazer(gmm.3.cont,gmm.3.cont.2,gmm.3.ntl,gmm.3.cw,gmm.3.orig,gmm.3.s,gmm.3.ns)


##############
###GED Data###
##############

#Remove Syria
main.grid.2 <- subset(main.grid, gwno!=652)

###Full
fel3.ged <- felm(nlights_mean~lagged_viol_tot+lag2ged_viol_tot+
                   laglogpop+lagnlights_mean+
                   excluded+p_polity2 +
                   logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid.2)

summary(fel3.ged)


####Full
##Keep only variables of interest
omit.full.ged <- na.omit(main.grid.2[,c("nlights_mean","lagged_viol_tot","lag2ged_viol_tot","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#omit.full$c_gid <- cluster(omit.full$gid)
#Convert data to pgmm formal
pgmm.full.ged <- pdata.frame(omit.full.ged, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.ged <- pgmm(nlights_mean~lagged_viol_tot+lag2ged_viol_tot+
                    laglogpop+lagnlights_mean+
                    excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                    lag(lagged_viol_tot, 2:7),
                  data = pgmm.full.ged, 
                  effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.ged)


################
###ACLED Data###
################

##Subset Africa and available temporal period (doesn't affect the results
main.grid.3 <- subset(main.grid, ((gwno>=400&gwno<=626)|gwno==651) )
main.grid.3 <- subset(main.grid.3, year>=1997)
###Full
fel3.acled <- felm(nlights_mean~lagacled_viol+lag2acled_viol+
                     laglogpop+lagnlights_mean+
                     excluded+p_polity2 +
                     logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid.3)

summary(fel3.acled)


####Full
##Keep only variables of interest
omit.full.acled <- na.omit(main.grid.3[,c("nlights_mean","lagacled_viol","lag2acled_viol","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#omit.full$c_gid <- cluster(omit.full$gid)
#Convert data to pgmm formal
pgmm.full.acled <- pdata.frame(omit.full.acled, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.acled <- pgmm(nlights_mean~lagacled_viol+lag2acled_viol+
                      laglogpop+lagnlights_mean+
                      excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                      lag(lagacled_viol, 2:7),
                    data = pgmm.full.acled, 
                    effect = "twoway", transformation= "d", model = "twostep")
summary(gmm.3.acled)

###Export to LaTex
stargazer(fel3.ged,gmm.3.ged,fel3.acled,gmm.3.acled)


#############################
###FEs and One-Way Effects###
#############################
###Year only FEs
fel3.a <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                       laglogpop+lagnlights_mean+
                       excluded+p_polity2 +
                       logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod
               |year|0|gid, data=main.grid)

summary(fel3.a)

####Gid only
fel3.b <- felm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                       laglogpop+lagnlights_mean+
                       excluded+p_polity2 +
                       logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod
               |gid|0|gid, data=main.grid)

summary(fel3.b)

#######One-way GMM effects
##Keep only variables of interest
omit.full <- na.omit(main.grid[,c("nlights_mean","lagall_perp_ml","lag2all_perp_ml","laglogpop","year","gid","lagnlights_mean","excluded","logwdi_expmil","logwdi_gdpcapcon2010","logross_oil_prod","logross_gas_prod","p_polity2")])
#omit.full$c_gid <- cluster(omit.full$gid)
#Convert data to pgmm formal
pgmm.full <- pdata.frame(omit.full, index=c("gid", "year"))

#set seed
set.seed(50)

###Run model
gmm.3.b <- pgmm(nlights_mean~lagall_perp_ml+lag2all_perp_ml+
                        laglogpop+lagnlights_mean+
                        excluded+p_polity2 +logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod+ logross_gas_prod| 
                        lag(lagall_perp_ml, 2:7),
                data = pgmm.full, 
                effect = "individual", transformation= "d", model = "twostep")
summary(gmm.3.b)

###Export to LaTex
stargazer(fel3.a,fel3.b, gmm.3.b)


########################
########################
###Summary Statistics###
########################
########################

#\textit{NTL}$_{it}$ 
c(min(main.grid$nlights_mean, na.rm=TRUE),median(main.grid$nlights_mean, na.rm=TRUE),mean(main.grid$nlights_mean, na.rm=TRUE),max(main.grid$nlights_mean, na.rm=TRUE),sd(main.grid$nlights_mean, na.rm=TRUE))

#\textit{NTL}$_{it-1}$ 
c(min(main.grid$lagnlights_mean, na.rm=TRUE),median(main.grid$lagnlights_mean, na.rm=TRUE),mean(main.grid$lagnlights_mean, na.rm=TRUE),max(main.grid$lagnlights_mean, na.rm=TRUE),sd(main.grid$lagnlights_mean, na.rm=TRUE))

#\textit{Atrocities}$_{it}$ 
c(min(main.grid$all_perp_ml, na.rm=TRUE),median(main.grid$all_perp_ml, na.rm=TRUE),mean(main.grid$all_perp_ml, na.rm=TRUE),max(main.grid$all_perp_ml, na.rm=TRUE),sd(main.grid$all_perp_ml, na.rm=TRUE))

#\textit{Atrocities}$_{it-1}$ 
c(min(main.grid$lagall_perp_ml, na.rm=TRUE),median(main.grid$lagall_perp_ml, na.rm=TRUE),mean(main.grid$lagall_perp_ml, na.rm=TRUE),max(main.grid$lagall_perp_ml, na.rm=TRUE),sd(main.grid$lagall_perp_ml, na.rm=TRUE))

#\textit{Atrocities}$_{it-2}$ 
c(min(main.grid$lag2all_perp_ml, na.rm=TRUE),median(main.grid$lag2all_perp_ml, na.rm=TRUE),mean(main.grid$lag2all_perp_ml, na.rm=TRUE),max(main.grid$lag2all_perp_ml, na.rm=TRUE),sd(main.grid$lag2all_perp_ml, na.rm=TRUE))

#\textit{Atrocities}$_{it-3}$ 
c(min(main.grid$lag3all_perp_ml, na.rm=TRUE),median(main.grid$lag3all_perp_ml, na.rm=TRUE),mean(main.grid$lag3all_perp_ml, na.rm=TRUE),max(main.grid$lag3all_perp_ml, na.rm=TRUE),sd(main.grid$lag3all_perp_ml, na.rm=TRUE))

#\textit{Atrocities}$_{it-4}$ 
c(min(lag(main.grid$lag3all_perp_ml,1), na.rm=TRUE),median(lag(main.grid$lag3all_perp_ml,1), na.rm=TRUE),mean(lag(main.grid$lag3all_perp_ml,1), na.rm=TRUE),max(lag(main.grid$lag3all_perp_ml,1), na.rm=TRUE),sd(lag(main.grid$lag3all_perp_ml,1), na.rm=TRUE))

#\textit{Atrocities (1 dic.)}$_{it-1}$ 
c(min(main.grid$lagall_perp_ml_one, na.rm=TRUE),median(main.grid$lagall_perp_ml_one, na.rm=TRUE),mean(main.grid$lagall_perp_ml_one, na.rm=TRUE),max(main.grid$lagall_perp_ml_one, na.rm=TRUE),sd(main.grid$lagall_perp_ml_one, na.rm=TRUE))

#\textit{Atrocities (10 dic.)}$_{it-1}$ 
c(min(main.grid$lagall_perp_ml_ten, na.rm=TRUE),median(main.grid$lagall_perp_ml_ten, na.rm=TRUE),mean(main.grid$lagall_perp_ml_ten, na.rm=TRUE),max(main.grid$lagall_perp_ml_ten, na.rm=TRUE),sd(main.grid$lagall_perp_ml_ten, na.rm=TRUE))

#\textit{Atrocities (20 dic.)}$_{it-1}$ 
c(min(main.grid$lagall_perp_ml_twenty, na.rm=TRUE),median(main.grid$lagall_perp_ml_twenty, na.rm=TRUE),mean(main.grid$lagall_perp_ml_twenty, na.rm=TRUE),max(main.grid$lagall_perp_ml_twenty, na.rm=TRUE),sd(main.grid$lagall_perp_ml_twenty, na.rm=TRUE))

#\textit{Atrocities (1 dic.)}$_{it-2}$ 
c(min(main.grid$lag2all_perp_ml_one, na.rm=TRUE),median(main.grid$lag2all_perp_ml_one, na.rm=TRUE),mean(main.grid$lag2all_perp_ml_one, na.rm=TRUE),max(main.grid$lag2all_perp_ml_one, na.rm=TRUE),sd(main.grid$lag2all_perp_ml_one, na.rm=TRUE))

#\textit{Atrocities (10 dic.)}$_{it-2}$ 
c(min(main.grid$lag2all_perp_ml_ten, na.rm=TRUE),median(main.grid$lag2all_perp_ml_ten, na.rm=TRUE),mean(main.grid$lag2all_perp_ml_ten, na.rm=TRUE),max(main.grid$lag2all_perp_ml_ten, na.rm=TRUE),sd(main.grid$lag2all_perp_ml_ten, na.rm=TRUE))

#\textit{Atrocities (20 dic.)}$_{it-2}$ 
c(min(main.grid$lag2all_perp_ml_twenty, na.rm=TRUE),median(main.grid$lag2all_perp_ml_twenty, na.rm=TRUE),mean(main.grid$lag2all_perp_ml_twenty, na.rm=TRUE),max(main.grid$lag2all_perp_ml_twenty, na.rm=TRUE),sd(main.grid$lag2all_perp_ml_twenty, na.rm=TRUE))

#\textit{Atrocities (non-ML)}$_{it-1}$ 
c(min(main.grid$lagall_perp_orig, na.rm=TRUE),median(main.grid$lagall_perp_orig, na.rm=TRUE),mean(main.grid$lagall_perp_orig, na.rm=TRUE),max(main.grid$lagall_perp_orig, na.rm=TRUE),sd(main.grid$lagall_perp_orig, na.rm=TRUE))

#\textit{Atrocities (non-ML)}$_{it-2}$ 
c(min(main.grid$lag2all_perp_orig, na.rm=TRUE),median(main.grid$lag2all_perp_orig, na.rm=TRUE),mean(main.grid$lag2all_perp_orig, na.rm=TRUE),max(main.grid$lag2all_perp_orig, na.rm=TRUE),sd(main.grid$lag2all_perp_orig, na.rm=TRUE))

#\textit{Atrocities (state)}$_{it-1}$ 
c(min(main.grid$lags_perp_ml, na.rm=TRUE),median(main.grid$lags_perp_ml, na.rm=TRUE),mean(main.grid$lags_perp_ml, na.rm=TRUE),max(main.grid$lags_perp_ml, na.rm=TRUE),sd(main.grid$lags_perp_ml, na.rm=TRUE))

#\textit{Atrocities (state)}$_{it-2}$
c(min(main.grid$lag2s_perp_ml, na.rm=TRUE),median(main.grid$lag2s_perp_ml, na.rm=TRUE),mean(main.grid$lag2s_perp_ml, na.rm=TRUE),max(main.grid$lag2s_perp_ml, na.rm=TRUE),sd(main.grid$lag2s_perp_ml, na.rm=TRUE))

#\textit{Atrocities (non-state)}$_{it-1}$ 
c(min(main.grid$lagns_perp_ml, na.rm=TRUE),median(main.grid$lagns_perp_ml, na.rm=TRUE),mean(main.grid$lagns_perp_ml, na.rm=TRUE),max(main.grid$lagns_perp_ml, na.rm=TRUE),sd(main.grid$lagns_perp_ml, na.rm=TRUE))

#\textit{Atrocities (non-state)}$_{it-2}$
c(min(main.grid$lag2ns_perp_ml, na.rm=TRUE),median(main.grid$lag2ns_perp_ml, na.rm=TRUE),mean(main.grid$lag2ns_perp_ml, na.rm=TRUE),max(main.grid$lag2ns_perp_ml, na.rm=TRUE),sd(main.grid$lag2ns_perp_ml, na.rm=TRUE))

#\textit{Atrocities (GED)}$_{it-1}$
c(min(main.grid.2$lagged_viol_tot, na.rm=TRUE),median(main.grid.2$lagged_viol_tot, na.rm=TRUE),mean(main.grid.2$lagged_viol_tot, na.rm=TRUE),max(main.grid.2$lagged_viol_tot, na.rm=TRUE),sd(main.grid.2$lagged_viol_tot, na.rm=TRUE))

#\textit{Atrocities (GED)}$_{it-1}$
c(min(main.grid.2$lag2ged_viol_tot, na.rm=TRUE),median(main.grid.2$lag2ged_viol_tot, na.rm=TRUE),mean(main.grid.2$lag2ged_viol_tot, na.rm=TRUE),max(main.grid.2$lag2ged_viol_tot, na.rm=TRUE),sd(main.grid.2$lag2ged_viol_tot, na.rm=TRUE))

#\textit{Atrocities (ACLED)}$_{it-1}$
c(min(main.grid$lagacled_viol, na.rm=TRUE),median(main.grid$lagacled_viol, na.rm=TRUE),mean(main.grid$lagacled_viol, na.rm=TRUE),max(main.grid$lagacled_viol, na.rm=TRUE),sd(main.grid$lagacled_viol, na.rm=TRUE))

#\textit{Atrocities (ACLED)}$_{it-1}$
c(min(main.grid$lag2acled_viol, na.rm=TRUE),median(main.grid$lag2acled_viol, na.rm=TRUE),mean(main.grid$lag2acled_viol, na.rm=TRUE),max(main.grid$lag2acled_viol, na.rm=TRUE),sd(main.grid$lag2acled_viol, na.rm=TRUE))

#\textit{Population}$_{it}$\textsuperscript{1} 
c(min(main.grid$laglogpop, na.rm=TRUE),median(main.grid$laglogpop, na.rm=TRUE),mean(main.grid$laglogpop, na.rm=TRUE),max(main.grid$laglogpop, na.rm=TRUE),sd(main.grid$laglogpop, na.rm=TRUE))

#\textit{Exclusion}$_{it}$
c(min(main.grid$excluded, na.rm=TRUE),median(main.grid$excluded, na.rm=TRUE),mean(main.grid$excluded, na.rm=TRUE),max(main.grid$excluded, na.rm=TRUE),sd(main.grid$excluded, na.rm=TRUE))

#\textit{Polity2}$_{jt}$ 
c(min(main.grid$p_polity2, na.rm=TRUE),median(main.grid$p_polity2, na.rm=TRUE),mean(main.grid$p_polity2, na.rm=TRUE),max(main.grid$p_polity2, na.rm=TRUE),sd(main.grid$p_polity2, na.rm=TRUE))

#\textit{Military expenditure}$_{jt}$\textsuperscript{1}
c(min(main.grid$logwdi_expmil, na.rm=TRUE),median(main.grid$logwdi_expmil, na.rm=TRUE),mean(main.grid$logwdi_expmil, na.rm=TRUE),max(main.grid$logwdi_expmil, na.rm=TRUE),sd(main.grid$logwdi_expmil, na.rm=TRUE))

#\textit{GDP PC}$_{jt}$\textsuperscript{1} 
c(min(main.grid$logwdi_gdpcapcon2010, na.rm=TRUE),median(main.grid$logwdi_gdpcapcon2010, na.rm=TRUE),mean(main.grid$logwdi_gdpcapcon2010, na.rm=TRUE),max(main.grid$logwdi_gdpcapcon2010, na.rm=TRUE),sd(main.grid$logwdi_gdpcapcon2010, na.rm=TRUE))

#\textit{Oil Production}$_{jt}$\textsuperscript{1}
c(min(main.grid$logross_oil_prod, na.rm=TRUE),median(main.grid$logross_oil_prod, na.rm=TRUE),mean(main.grid$logross_oil_prod, na.rm=TRUE),max(main.grid$logross_oil_prod, na.rm=TRUE),sd(main.grid$logross_oil_prod, na.rm=TRUE))

#\textit{Gas Production}$_{jt}$\textsuperscript{1}
c(min(main.grid$logross_gas_prod, na.rm=TRUE),median(main.grid$logross_gas_prod, na.rm=TRUE),mean(main.grid$logross_gas_prod, na.rm=TRUE),max(main.grid$logross_gas_prod, na.rm=TRUE),sd(main.grid$logross_gas_prod, na.rm=TRUE))

#\textit{Drought}$_{it}$ 
c(min(main.grid$droughtyr_spi, na.rm=TRUE),median(main.grid$droughtyr_spi, na.rm=TRUE),mean(main.grid$droughtyr_spi, na.rm=TRUE),max(main.grid$droughtyr_spi, na.rm=TRUE),sd(main.grid$droughtyr_spi, na.rm=TRUE))

#\textit{Prec.}$_{it}$\textsuperscript{1} 
c(min(main.grid$logprec_gpcp, na.rm=TRUE),median(main.grid$logprec_gpcp, na.rm=TRUE),mean(main.grid$logprec_gpcp, na.rm=TRUE),max(main.grid$logprec_gpcp, na.rm=TRUE),sd(main.grid$logprec_gpcp, na.rm=TRUE))

#\textit{Temp.}$_{it}$ 
c(min(main.grid$temp, na.rm=TRUE),median(main.grid$temp, na.rm=TRUE),mean(main.grid$temp, na.rm=TRUE),max(main.grid$temp, na.rm=TRUE),sd(main.grid$temp, na.rm=TRUE))

#\textit{UCDP CW}$_{jt}$ 
c(min(main.grid$ucdp_type3_bin, na.rm=TRUE),median(main.grid$ucdp_type3_bin, na.rm=TRUE),mean(main.grid$ucdp_type3_bin, na.rm=TRUE),max(main.grid$ucdp_type3_bin, na.rm=TRUE),sd(main.grid$ucdp_type3_bin, na.rm=TRUE))

#\textit{UCDP CW}$_{jt-1}$ 
c(min(main.grid$lagucdp_type3_bin, na.rm=TRUE),median(main.grid$lagucdp_type3_bin, na.rm=TRUE),mean(main.grid$lagucdp_type3_bin, na.rm=TRUE),max(main.grid$lagucdp_type3_bin, na.rm=TRUE),sd(main.grid$lagucdp_type3_bin, na.rm=TRUE))

#\textit{Cap. distance}$_{it}$\textsuperscript{1}  
c(min(main.grid$logcapdist, na.rm=TRUE),median(main.grid$logcapdist, na.rm=TRUE),mean(main.grid$logcapdist, na.rm=TRUE),max(main.grid$logcapdist, na.rm=TRUE),sd(main.grid$logcapdist, na.rm=TRUE))

#\textit{Arms imp.}$_{jt}$\textsuperscript{1} 
c(min(main.grid$logwdi_armimp, na.rm=TRUE),median(main.grid$logwdi_armimp, na.rm=TRUE),mean(main.grid$logwdi_armimp, na.rm=TRUE),max(main.grid$logwdi_armimp, na.rm=TRUE),sd(main.grid$logwdi_armimp, na.rm=TRUE))

#\textit{Life exp.}$_{jt}$ 
c(min(main.grid$wdi_lifexp, na.rm=TRUE),median(main.grid$wdi_lifexp, na.rm=TRUE),mean(main.grid$wdi_lifexp, na.rm=TRUE),max(main.grid$wdi_lifexp, na.rm=TRUE),sd(main.grid$wdi_lifexp, na.rm=TRUE))

#\textit{CIRI}$_{jt}$
c(min(main.grid$ciri_empinx_new, na.rm=TRUE),median(main.grid$ciri_empinx_new, na.rm=TRUE),mean(main.grid$ciri_empinx_new, na.rm=TRUE),max(main.grid$ciri_empinx_new, na.rm=TRUE),sd(main.grid$ciri_empinx_new, na.rm=TRUE))




fel3.cal <- felm(nlights_calib_mean~lagall_perp_ml+lag2all_perp_ml+
               laglogpop+lag(nlights_calib_mean, 1)+
               excluded+p_polity2 +
               logwdi_expmil+logwdi_gdpcapcon2010+logross_oil_prod + logross_gas_prod|year+gid|0|gid, data=main.grid)

summary(fel3.cal)